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We show that in anhamionic one-dimensional crystal lattices pairing of electrons or holes in a 
localized bisolectron state is possible due to coupling between the charges and the lattice defor- 
mation that can overcompensate the Coulomb repulsion. Such localized soliton-like states appear 
as traveling ground bound singlet states of two extra electrons in the potential well created by 
the local lattice deformation. We also find the first excited localized state of two electrons in a 
soliton-Uke lattice deformation potential well given by a triplet state of two electrons. The results 
of the analytical study of interacting electrons in a lattice with cubic anharmonicity are compared 
with the numerical simulations of two electrons in an anharmonic lattice taking into account of the 
(local) Hubbard electron-electron repulsion. We qualitative agreement between both approaches for 
a broad interval of parameter values. For illustration we give expressions for the bisolectron binding 
energy with parameter values that are typical for biological macromolecules. We estimate critical 
values of Coulomb repulsion where the bisolectron becomes unbound. 

PACS numbers: 71.10.Li, 63.20.Ry, 71.38-k, 71.38.Mx 



I. INTRODUCTION 

A broad class of low-dimensional systems such as poly- 
diacetylene [l|-lj] , conducting platinum chain compounds 
and conducting polymers [5|, salts of transition metals 
(PbSe, PbTe, PbS) [l-i], superconducting cuprates [H- 
|l5| which find numerous applications in microelectronics 
and nanotechnologies, or play important role in li ving 
systems (polypeptide macromolecules, DNA, etc.) |16l - 
|23|, manifest nonlinear effects. These effects arise from 
the electron-lattice interaction or from the (nonlinear) 
lattice anharmonicity and, eventually, from both items 
coupled together. As a consequence there is the forma- 
tion of soliton-like states, that have been called electro- 
solitons or solectrons in different albeit related contexts 
[m |24| . Both generalize the polaron concept |25l - [3l| . It 
has been shown that anharmonic lattices can support the 
formation of two-component solitoi is | 32h 34| and pairing 
of two excess electrons or holes |35l - [39| may be enhanced 
bringing strong correlation both in momentum space and 
in real space. 

The problem for pairing of charged particles is the 
Coulomb repulsion. In normal metallic superconductors, 
the attractive interaction needed to form a bound state 
of electrons (Cooper pair) is due to phonon excitations in 
the harmonic lattice. Lattice polarization overcomes the 
Coulomb repulsion so that for states near to the Fermi 
surface, where the repulsion is rather weak, pairing be- 
comes possible (in momemtum space only with complete 
delocalization in real space). This situation changes in 
the strong coupling case where the two-electron bound 
state is localized in space on one lattice site (small bipo- 



laron) ^,^. 

For one-dimensional electron-lattice systems with cu- 
bic or quartic anharmonicity, we have shown |38l . |39[ that 
excess electrons with opposite spins may form a bisolec- 
tron, which is a localized bound state of the paired elec- 
trons in a soliton-like lattice deformation. When one ex- 
cess electron is localized on several lattice sites forming a 
Solectron state, that is a bound electron-soliton state, the 
Coulomb repulsion between two electrons in a bisolectron 
is not as strong as in a small bipolaron, but is not negli- 
gible as in the case of completely dclocalized electrons. 

In contrast to two- and three-dimensional systems, the 
Coulomb interaction can lead to divergent potential en- 
ergy in one-dimensional systems if both electrons are 
found with finite density at the same point in space. 
Nevertheless, this singularity can be removed by using 
the Hubbard model of a discrete lattice [40, |41|. We 
consider in Sec. II the case when the electron-lattice in- 
teraction is moderately strong, so that an adiabatic ap- 
proximation is valid which leads to the soliton-like solu- 
tions of the system of coupled nonlinear differential equa- 
tions in the continuum approximation. Then we study 
first the case when the Coulomb interaction is relatively 
weak thus modifying weakly the two-electron bisolectron 
wave function (Sec. Ill) and, subsequently, the rather 
strong Coulomb repulsion case for which we find the 
two-electron wave function within a variational approach 
(Sec. IV). Solutions for the ground state and the first ex- 
cited state are considered in Sec. V. We compare these 
analytical results with results of numerical simulations 
for two excess electrons in the lattice in Sec. VI. The 
lattice interaction between the nearest neighbors is de- 



scribed by the Morse potential that gives also a cubic 
anharmonicity. The electrons are considered within a 
Hubbard model. Finally, in Section VII we draw some 
conclusions, discussing the applicability of the continuum 
model and the stability of the bisolectron state. 



II. DYNAMICS OF THE COUPLED 
ELECTRON-LATTICE SYSTEM 

We consider an infinitely long, one-dimensional (Id) 
crystal lattice, with units all of equal mass M and equi- 
librium lattice spacing a, where two free, excess electrons 
are added. The Hamiltonian of such a system can be rep- 
resented in the form: 



H — Hf,\ + H] 
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Let Eq denote the on-site electron energy, let J denote 
the electron exchange interaction energy, and let cjj ^ and 
Cn.s be the creation and annihilation operators of an elec- 
tron with the spin projection s =t, i at the site n. Then 
the electron Hamiltonian in ([1]) has the explicit form 
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With the account of only longitudinal displacements of 
atoms from their equilibrium positions, the lattice part 
of the Hamiltonian in ^ has the form 
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where u„ is an operator of the n-th atom displacement, 
Pn is the operator of the canonically conjugated momen- 
tum, [/lat is the operator of the lattice potential energy, 
whose properties will be defined below. 

We consider the case when the dependence of the 
on-site electron energy on atom displacements is much 
stronger than that of the exchange interaction energy, so 
that the electron-lattice interaction part of the Hamilto- 
nian ([!]) has the form 
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where x is the electron- lattice interaction constant. 

The Coulomb repulsion between the electrons is given 
by the Hamiltonian 
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(5) 



where C/„„i is the matrix element corresponding to the 
Coulomb interaction of electrons on sites n and m. 

In the Born-Oppcnheimcr approximation we can fac- 
torize the state vector of the system into two parts. 



\^{t))^\^,,{t))\^ut{t)). 



(6) 



Here the state vector of the lattice has the form of the 
product of the operator of coherent displacements of the 
atoms and the vacuum state of the lattice 

|*lat(t)) = 5'(t)|0)lat, 

S{t) == CXp S -tE Wn{t)Pn ~Pn{t)Un] \ , (7) 



where Un{t), Pn{t) are, respectively, the mean values of 
the displacements of atoms from their equilibrium posi- 
tions and their canonically conjugated momenta in the 
state dH). 

In the case of two excess electrons with spins si, S2 in 
the lattice, the electron state- vector has the form 

l*ci(i)) = Yl *("i'^2,si,s2;t)c].^^,^cj,^_,j0)ci,(8) 

where the function ^(ni, n2, Si, S2; i) satisfies the nor- 
malization condition 



^ |*(ni,n2,Si,S2;t)P = l- 
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In the absence of a magnetic field the electron wave 
function of the system can be represented by the product 
of the spatial and the spin wave functions. The antisym- 
metry requirement of the two-electron wave function may 
be fulfilled by the symmetry of the spatial wave function 
and the antisymmetry of the spin function (singlet state) 
or by the antisymmetry of the coordinate function and 
the symmetry of the spin wave function (triplet state). 
Using the state vector ([SI)-® we can calculate the Hamil- 
tonian functional iJ, corresponding to the Hamiltonian 
operator ([1]): 



H^{^{t)\H\^{t)). 



(10) 



Minimizing the functional (fTUI) with respect to electron 
and phonon variables, we derive the system of two cou- 
pled equations for the two-electron wave function and 
the atom displacements. In the continuum approxima- 
tion n ^>- X = na this system of equations has the form: 



_ d-9{xi,X2,t) 
' dt 

fdu{x,t) 



2m \dxl 
du{x,t) 



dxl 



+t/coui(a;i,a;2)*(a;i,X2,t), 



dt^ ^' ap2 dx^ 



^{xi,X2,t) 

*(.Ti,X2,i) 

(11) 



' M^ [dxT / ^^2|^(a;i,a;2,i)PUi=:i 

-^ / dxi\-^{xi,X2,t)\'^\x2=x] ■ 



(12) 



Here, 14c — a^JwjM is the sound velocity of the lat- 
tice, w is the lattice elasticity coefficient, \['(xi,X2,i) is 
the two-electron spatial wave function with xi = nia^ 
i = 1,2, u{x,t) is the function of the lattice displace- 
ment, and C/coui accounts for the Coulomb repulsion be- 
tween the two electrons. Note that to study the case of 
strong electron localization we have to include into equa- 
tion (fT^ an additional term ex d'^u/dx^ which takes into 
account the dispersion of the lattice and can lead to the 
supersonic solutions, see, e.g., |42| . 

Introducing the function of the lattice deformation 
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we can write the potential energy of the system appearing 
in the right hand side of the Schrodinger equation (jllI) 
in the form 

C^tot = C/soi(a;i,X2,i) + C/coui(a;i,X2), (14) 

Uso\ixi,X2,t) = -xa [p(a;i,i) + pix2,t)] , (15) 

where the function Usoi{xi,X2,t) plays the role of an ef- 
fective potential created by the deformation of the lattice 
(index 'sol' indicates that this potential is a soliton-like, 
as it will be shown below). 

For the case of the lattice with cubic anharmonicity 
the potential energy of the lattice, measured in units of 
J, is 
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where 7 is the anharmonicity parameter whose value is 
positive, as we wish to focus on rather strong compres- 
sions for which the repulsion part of the potential is what 
really matters. 

The solution of Eqs. ([TT|l - p2|) without account of the 
Coulomb repulsion, i.e., when in Eq. PT|) t/coui is omit- 
ted, was found previously [3M| for the singlet two-electron 
problem. In particular, it was found that the lattice de- 
formation is given by the expression 



p{0 = Posech^KO, ^ = {x-Vt)/a, (17) 



For arbitrary values V of the velocity, po{V) can be found 
numerically. 

The solution for the spin-singlet two-electron wave 
function is 

*(a;i, .T2, t) = -= [^i{xi,t)^2ix2,t) + ^2{xi,t)^iix2, t)] 

(21) 
The single-electron wave functions \['i(a;,i) = 'i>2{x,t) = 
^(a;,i) are normalized to 1, according to Ec^. ([9]). The 
modulus [^^(a;,^)! = $(^) is found to be of the soliton- 
type form 
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with s^ = V'^ /V^^. The phase factor contains the energy 
and will no be given here. Thus, the parameter k in Eqs. 
(fT7)) and ([m determines the width of electron localiza- 
tion in the bisolectron state (|22|) . formed due to binding 
of the two electrons with the lattice deformation in the 
soliton state (|17|) without account of the Coulomb repul- 
sion. As shown in [38| , the singlet bisolectron is always a 
bound state compared with two independent solectrons 
if the Coulomb repulsion is neglected. The following re- 
lation for the binding energy has been found: 
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(23) 
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Here, E^^°^'> and Pq"'' are the energy and maximum lat- 
tice deformation for the solectron state where one elec- 
tron is bound to the lattice soliton. The corresponding 
values are given in Ref. [38|. For instance, ty pica l val- 
ues for macromolecules like polypeptides (see [131) ^-''s: 
X = 35 - 62 pN, a == 5.4 • 10"^" m, u; === 39 - 58 N/m, 
M = 5.7-10-25 kg, Kc = (3.6-4.5)-103m/s. The trans- 
fer energy value for the amid-I excitation is Jexc = 0.001 
eV, and the the electron transfer energy is of the or- 
der J « 0.1 — 0.5 eV. For these parameter values, the 
bisolectron binding energy ((23|) in macromolecules can 
be estimated in the region ijt^™'*) « 0.05 — 0.5 eV, if the 
Coulomb repulsion is neglected. 
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III. LIMIT CASE OF WEAK COULOMB 
REPULSION 



The maximum value of the lattice deformation po de- 
pends on the soliton velocity V and lattice anharmonic- 
ity. For the stationary case, 1^ = 0, it is given approxi- 
mately by the analytical expression 
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If the bisolectron state is extended over few lattice 
sites, i.e., if the width of the bisolectron l^^^'^^a ~ 27ra/K 
is bigger than a, the Coulomb repulsion is relatively weak 
as comparing with the binding energy of the bisolectron 
(P^ . Therefore, the wave function (P^ can be general- 
ized as (see [3^, [Sy, [33] ) 



^i.AO 




(24) 



xsech[K(^ =F V2)] V 1 - s^ + -fPosech^[K{^ T 1/2)] . 



Here the parameter I accounts for the Coulomb repulsion 
in Eq.(|ll|). so that approximately 
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and determines the distance between the centre of mass 
(c.o.m.) coordinates of one-electron wave functions. Here 
e denotes the dielectric constant of the lattice. Note that 
here and below we use SI units, and e = Creieo contains 
the relative dielectric constant erci of the medium. 

The parameter / can be determined from the condition 
of the minimum of energy. In the weak Coulomb repul- 
sion case where 2tt/k ^ a, the total energy of the system 
is 
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Minimizing this expression with respect to /, we get the 
equilibrium distance between the maxima of one-electron 
functions: 
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Then (j27p can be approximated by 
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We remind here that this approach strictly speaking, is 
valid for relatively weak Coulomb repulsion, so that the 
condition Iq < /('"'"^ = 27r/K is fulfilled. Thus we find 
a lattice soliton binding two electrons together [Sgl in 
the lattice deformation potential well p3|) . The bisolec- 
tron solution is stable when its binding energy is larger 
than the Coulomb repulsion that means with Eq. ((23l) 
Jg'^/2 > C/coui, where g = x^/C^Jw) = (0.9 - 3.9) for 
polypeptides [13 ■ Therefore we conclude that the critical 
value of the electron-lattice coupling constant increases 
with increasing Coulomb repulsion that agrees with the 
results of the numerical simulations given in Ref. [5l| . 
Note that the Coulomb interaction is screened by the 
relative dielectric constant e^c\ of the order of 10. For the 
parameter values given above, we find that Iq is small 
compared with 2tt/k so that the weak coupling condition 
holds. The estimation for the binding energy E^^^^^ gives 
a value of about 0.05 eV indicating a stable bisolectron 
solution. 

In Section VI, FigHl we will show the charge density 
distribution, q{£,) = $?(0 + *2(0 (sec Eq. ^^) for three 
different values of Iq (sec (P7|) or ^^), in other words. 



for three different relations between the binding energy 
of the bisolectron and Coulomb repulsion ([25|) . We will 
see there that the density profile has one maximum at 
relatively small values of Iq or two maxima at large values 

of /o- 

These results, as will be shown in Section VI, with high 
degree of accuracy explain the results of numerical sim- 
ulations of two electrons in a discrete Morse lattice with 
account of Hubbard interaction between the electrons in 
a broad interval of the strength of interaction. 



IV. SEPARATION OF THE TWO-ELECTRON 
CENTER OF MASS MOTION 

Here we consider the case of arbitrary Coulomb re- 
pulsion, when the Coulomb term in Eq. pip can be 
represented as 



C^Coul = 



ATTe\xi - X2\ 



(30) 



Let us introduce now the c.o.m. coordinate X ~ {xi + 
x-i)l2 — Vt, the running wave coordinate of the soli- 
ton moving with velocity V , and the relative coordinate 
X = a:2 — xi . The interaction of electrons with the lattice 
is described by the effective deformational potential (|15p . 
that, according to the results found in the previous Sec- 
tion, is created by the soliton and may bind the two elec- 
trons together. We are assuming now that the solution 
for the soliton potential C/soi(a;i, X2, t) can be represented 
in the form: 

Usol{xi,X2,t) « U!^"{{xi,X2,t) 



-Xa 



p('\xi,t)+p('»ix2,t) 
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(31) 



where superscript (0) indicates solutions without ac- 
count of the Coulomb interaction, and where we expand 
Urci{x,t;X) to second order with respect to x. We use 
the approximation 
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Since the potential (|3ip splits into two parts we may sep- 
arate the Hamiltonian functional 
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Here we have used only the first non- vanishing term of the 
expansion of the potential around its minimum. Higher 
orders can be treated by perturbation theory. In general. 




Figure 1: Deformation potential (|35p for 7 = 1.5, S = 0.002 
po ~ 0.007 and k — 1.2, (black solid line) and a Gaussian (|36( 
fitted to the region around the minimum (red dashed line) 
Potentials are given in units 2xa,po- 



the relative potential well in the co-moving coordinates 
Uiei{x-,t; X) depends not only on the relative coordinate 
X but also on the center of mass position X so that the 
relative motion and the c.o.m. motion become coupled. 

Our aim was to separate the problem into a c.o.m. 
problem identical to the previous task without Coulomb 
repulsion and a parabolic problem for the c.o.m. motion. 
This way we have first to solve the soliton problem for 
the common motion to find the function p^"^ {X, t) and 
Wo, and then to solve the parabolic problem including the 
Coulomb term (pO|) . 

With the ansatz '^{xi,X2,t) = $(X)0(a;), E = 
Ec.o.m. + E^ci, we get the equation for the center of mass 
motion 



'ima^ dS,^ 
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HaO = E,,o.n.MaO- (34) 



where ^ — X/a. 

Here the soliton deformation p^^\C) is determined in 
Eq. p7|) . and, therefore the soliton part of the deforma- 
tion potential ([5T|) is given by expression 



C^soi(0 = -2xa/0osech^«). 



(35) 



We can approximate expression p5[) by the following one 



UcsiO = -2xapoe 
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(36) 



that coincides with expression ([35|) up to the order ^^. 
For illustration in Fig. [1] we show the deformation po- 
tentials for the case 7 = 1.5 and other parameters corre- 
sponding to their values in poly-peptide macromolecules. 
For the given parameter values the deformation poten- 
tials given by expressions ([35]) and (|36p show negligi- 
ble difference, especially near the minimum, which de- 
termines the potential C/rci ((32|) of relative motion. The 
parabolic approximation ([5^ is possible only for x < 
a/n. For larger distances, the deformation potential be- 
comes weak. 



V. WAVE FUNCTION FOR THE 

TWO-ELECTRON RELATIVE MOTION 

INCLUDING COULOMB REPULSION 

Let us now solve the Schrodinger equation for the rel- 
ative motion of electrons: 






47re a; 



(j)ix) = E,,i(bix) (37) 



with the symmetry condition for fermions that in the 
spin singlet state 0(— x) = (j){x) (symmetric orbital) and 
in the spin triplet state (l){—x) = —4>{x) (antisymmetric 
orbital). 



The minimum of the potential follows from p7p at x 
xq, where 
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(38) 



that may be considered as a characteristic length scale. 
In a classical treatment, it gives the distance between the 
two electrons at zero temperature. 

In the absence of Coulomb repulsion, the solutions for 
the harmonic potential are Hermitean functions. The 
ground state is given by the symmetric orbital that has 
no nodes, but has a finite density at x = that leads 
to a divergent energy if Coulomb repulsion is considered. 
Therefore, with account of Coulomb repulsion we have in 
the Id case the additional condition 0(0) = for strongly 
localized electrons, i.e., only solutions with a node at a; = 
are possible. To consider the behavior near a; = 0, we 
have (j){x) oc x + e^m/ {2h'^e)x^ . Therefore, both electrons 
are kept apart by a distance of the order a;o given by the 
value (p8|) . Comparing this expression with the result 
([27]) . it should be mentioned that now the deformation 
of the lattice remains fixed when the distance between 
the electrons changes. 

A rigorous solution can be given in the free case wq = 
using hypergeometric functions 



x) = const xexp —i{mEi-c\/h ) ' x 



(39) 
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This solution approximates the full solution with the 
harmonic potential (or another potential that is finite 
near a; = 0). Including the harmonic potential, the 
full solutions can be found numerically. As a peculiar- 
ity of the one-dimensional Coulomb problem, symmet- 
ric solutions 4>{—x) = 4'{x) and antisymmetric solutions 
(t>{—x) = —4>{x) become degenerated. This is an artifact 
that arises performing the continuum limit for strongly 
localized solutions and is removed if we consider the dis- 
crete lattice and take into account the higher dimension- 
ality of real systems (even so-called one-dimensional sys- 
tems are in fact confined three-dimensional systems). 



In [38| we discuss the solution (f>{x) of the Schrodinger 
equation for the relative motion using a variational ap- 
proach. Taking into account the properties we already 
discussed, we consider the normalized function 
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that is the Hermitean function for the first excited state. 
We can improve this ansatz using, e.g., higher Hermitean 
functions or other zero-node functions like a;sech(/3ix). 
The parameter j3 is now a free variational quantity. In the 
variational problem the energy is the sum of the kinetic 
energy, harmonic potential and the Coulomb part. For 
the class of functions (1^0)1 we have 



E,[p)=Z^P+^-^ + J-—P^'\ (41) 
m 8p V TT 47re 



The condition for the minimum gives the minimal energy 
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In the limit where the Coulomb term can be neglected, 
we find fi\ = muJo / {2y/2h) . In the opposite case 

where the Coulomb part is strong, we have /3i = 

2/3 



2 8"?" 



— 1/2 

, so that 13-^ is proportional to xq, 



Eq. (|38|). As an interpolating formula we can take 
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A more sophisticated ansatz for the variational treat- 
ment would be the following class of functions 






^ g-/3(x-o' _ e"'^("^+') 



(44) 



In the limit I — ?> the Hermitean function (|40p is recov- 
ered. The corresponding energy functional £'|'^(/3,/) = 



ffr, 



for the Hamiltonian of the relative motion 



can also be given in analytical form that is very tedious 
and will be discarded here. 

With respect to the symmetric solution, we can take 
\(t>i{x] /3)| from Eq. ((40|) that gives the same energy value 
as the antisymmetric one. A symmetric wave function 
0|(a;; /3, /) that has no nodes can be obtained as sym- 
metric superposition similar to Eq. (j44p . However, the 
energy is diverging because 0|(O;/?, ?) 7^ 0. If the ar- 
tificial singularity of the Coulomb repulsion is removed 
considering the distribution of the electron wave function 
in the space as given by the atomic orbits, the symmetric 
solution may become favorable. For instance, we can re- 
place the Coulomb repulsion by e? /{A'K€\/x'^ -\- a^) where 



Figure 2: Potential energy of the relative motion (in Rydberg 
units) and wave function for the relative motion of the elec- 
trons from the numerical solution of the Schrodinger equation 
(|37p . for different parameter values of ujq — mojo/Ryd. 



a is of the order of the Bohr radius. Then, the degenera- 
tion between the symmetric and antisymmetric solution 
will be removed so that the symmetric solution (spin sin- 
glet) becomes energetically favorable. We will discuss 
this point below in context with the Hubbard model. 

Instead of the variational approach, the Schrodinger 
equation ([37]) for the relative motion can be solved nu- 
merically. Potentials for different parameter values of ujq 
(in Rydberg units, h = e^ /{^ire) = 2m = 1) are shown 
in Fig. [21 The corresponding normalized wave functions 
are also shown in Fig. [2] for a; > 0. For a; < we have the 
symmetric solution (f){—x) = 4>{x) and the antisymmetric 
solution (j){—x) = —4>{x) that are degenerated. 

We see that the variational ansatz (|40p can be used 
and coincides well with the exact solution. However, we 
point out that the harmonic approximation for the de- 
formation potential is justified only near the minimum of 
the potential, i.e. x < a/n. 

To consider also larger value of the relative distance 
X, we can approxiinate the interaction by a potential of 
solitonic form, instead using the harmonic potential, 

e^ \ 

usech (kx) + rrr ) (j^ix) = £'rcl'/'(S). 

47rea|x| ) 

(45) 

The relative distance x = xj a is given in units of a. The 
solution of this equation indicates the range of parame- 
ter values where a bisolectron solution exists. With the 
parameter values given for polypeptides [l7| , we estimate 
the values h^ /ma? « 0.2 eV, e^ /Airea k. 0.5 eV. A bound 
state exists if u exceeds a critical value Uc- As example. 



dx"^ 



we find Uc = 0.4742 cV for k = 0.3, Uc = 0.1368 cV for 
R = 0.1. Wc can estimate R ~ k/\'2,u w X'^Po- If "^^ 
consider the value for _b(*''"'^) given at the end of Sec. |Tl] 
to characterize the parameter range of u, a bound bisolec- 
tron state is stable for values of R below the correspond- 
ing critical values that are in the region of parameter val- 
ues relevant for polypeptides. Note that our approach, 
which considers the relative motion of the two electrons 
separately from the center-of-mass motion, works also 
in the Hubbard models for the bisolectrons considered 
in p3l - |45| . There, the potential strength —2axPo is re- 
placed by —Ue that is the minimum of the polarization 
potential. 



VI. COMPARISON WITH NUMERICAL 
RESULTS FOR THE HUBBARD MODEL 



Our treatment of bisolectron states in one-dimensional 
anharmonic lattices contains two basic approximations. 
To obtain analytical expressions, we considered the dis- 
crete lattice in continuum approximation. Because the 
extension of the bisolectron state is not very large, the 
discrete structure of the lattice may become of relevance. 
In particular, it is not simple to relate the parameter 
values of the continuum model with the properties that 
describe the lattice. Thus, the Coulomb repulsion is over- 
estimated a short distances in the continuum model, and 
the singularity that occurs at zero distance between the 
two electrons is an artifact. This singularity is removed 
when the Hubbard model is introduced. The other ap- 
proximations that can be controlled are the expansions 
in deriving the nonlinear potential and the separation of 
the center of mass motion from the relative motion. For 
the parameter values considered here, these approxima- 
tions seem to be reasonable and can be improved adding 
further terms of the expansions. 

We will compare now our analytical results with sim- 
ulations based on the Hubbard Hamiltonian, coupled to 
a discrete lattice with Morse interaction [4J, |4a] (for re- 
lated work see |47H51| ). The Hamiltonian we use for the 
simulations describes a Id periodic lattice of molecules 
in which two excess electrons have been injected. For 
estimation we give some typical parameter values of the 
Morse potential 



U 



Morse 



(r) =i:)[(l-e--^('-''))2-l] 



(46) 



that are adapted to describe biomolecules: the lattice 
spacing a ~ 1 — 5 A, the stiffness Ba ~ 1, the depth 
D ~ 0.1 — 0.5 eV. The period of oscillations in the Morse 
well results as l/r^MorsG — 0.1 — 0.5 ps. 

Noteworthy is that for the Toda potential [54I, with 
about the same repulsive component, we know the exact 
analytical solutions that, for the Morse potential (|46)) . 
can be used to quite a satisfactory level of approxima- 
tion [2J, |5j, |5J] . This potential energy of a lattice with 
Morse interaction can be expanded in a power series with 



respect to the displacements of the atoms from equilib- 
rium positions. The terms beyond the quadratic one de- 
scribe the anharmonicity of the lattice. Here, we are 
interested in the cubic term that can be compared with 
our analytical approach. 

The electron part is given by the Id Hubbard Hamil- 
tonian, cf. Eqs. (HI), (0), 



Hcl + -ffcoui 



+ UY.^lAi^ni£nf (47) 



Instead of the Coulomb interaction also between dis- 
tant lattice sites, only the single-site Hubbard repulsion 
Un.m = USn,m is Considered. The quantities Jn,n+i de- 
note the transfer matrix elements whose value are de- 
termined by an overlap integral being responsible for 
the nearest-neighbor transport of the electron along the 
chain. These transfer matrix elements are supposed here 
to depend on the relative distance between two consecu- 
tive units in the following exponential fashion |55| 



Jn,7i+i = Jo exp[-?7 (u„+i - u„)] ■ 



(48) 



The quantity 77 regulates how strong the Jn,n+i are influ- 
enced by the relative displacement of lattice units, or in 
other words it determines the coupling strength between 
the electron and the lattice system. The numerical meth- 
ods to solve the Hubbard model system coupled to the 
Morse lattice are given in Refs. |46l.[5q|. 

Compared with the continuum model used in this 
work, the Coulomb repulsion between the electrons is de- 
scribed in the Hubbard model only in quite ad-hoc albeit 
appropriate approximation. We use the comparison with 
the Hubbard model to control the approximations per- 
formed in our approach. In particular, the transition to 
the continuum limit is possible when the relevant prop- 
erties vary smoothly in the lattice. With respect to the 
Coulomb interaction this means that the characteristic 
length scale given here by the width of the bisolectron 
1/k is large compared with the lattice parameter. In the 
opposite case, the simple Coulomb interaction that is sin- 
gular at a; = cannot be used. A Hubbard Hamiltonian, 
in general with two-center Coulomb repulsion terms, may 
become more appropriate. 

The numerical results for electron pairs with single-site 
Hubbard repulsion are shown in Fig. [3] for the electron 
pair density (blue solid line) and the particle velocity dis- 
tribution (red dashed line). The parameter values used 
in the simulations were ?/ — 2.5a, Jq = 0.041?. The adi- 
abaticity parameter r = Jo/(^^Morsc) was fixed at the 
value T = 20. 

We discuss the parameter values that are in the typ- 
ical range to describe biomolecules more in detail. The 
Hubbard parameter was also measured in units of the 
oscillation energy of the Morse potential using the no- 
tation U ~ U/IiQmoi-sc- We consider in the numeri- 
cal simulations the values U — 20, 60, 100. We used 
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Figure 3: Numerical estimates for the electron density (blue 
solid line) and the velocity distribution (red dashed line) of 
Solectron pairs with Hubbard repulsion at (7 = 20, 60, 100 on 
the Morse lattice. 



Figure 4: Estimates of the bisolectron charge density profile 
9(0 1 "^s. the lattice coordinate ^ using the analytical result 
((221) with 7 = 1.5, K = 0.35, and h = 2, 6, 10. 



S^Morse = 0.0021?. With values of the depth of the 
well in the range of D ~ 0.1 — 0.5 eV we see that the 
lowest f7 = 20 would correspond to C/ = 0.004 - 0.02 
eV and the upper value U = 100 would correspond to 
\J — 0.02 — 0.1 eV. Coulomb repulsion energies in this 
range is what one expects from physical estimates for 
electron pairs in a medium with a relative dielectric con- 
stant around Croi — 10. Summarizing we may conclude 
that Hubbard parameters in the range C7 = 20 — 100 
which we considered here correspond to Coulomb repul- 



sion energies in the range U ~ 0.01 — 0.05 eV what seems 
to make sense. 

We see, that depending on the repulsion there can be 
one or two peaks in the both characteristics, the probabil- 
ity distribution and the velocity distribution. The proba- 
bility distribution gives the probability to find an electron 
at the lattice point under consideration, independent on 
the position of the other electron. It measures the elec- 
tron density. The velocity distribution proofs that there 
is a collective flow as expected for a solitonic excitation. 
We emphasize that the simulation also shows the stabil- 



ity of the bisolectron that stays bound during the whole 
run of shnulation. 

In Fig. m we show the charge density function within 
our analytical model developed in Section III, for various 
values of the Coulomb repulsion, which determines the 
distance between the maxima of one-electron functions. 
Note that the charge density is determined by (|25|) with 
I = Iq ((77)) . We set the ratios of Iq in the three figures 
the same as the ratios of U in Fig. [3l 

We cannot expect a one-to-one correspondence be- 
tween the the numerical and analytical results. They 
are obtained in different models to treat the anharmonic 
lattice and the Coulomb repulsion. However, we see a 
qualitative agreement in both approaches. In particu- 
lar, we see that electrons are localized in the bisolectron 
state, the profile of which depends on the strength of the 
Coulomb repulsion. The tendency that one maximum 
splits into two maxima with increasing Coulomb repul- 
sion is clearly seen. 

Analyzing the results shown in Figs. [31 [5] we see that 
increasing the Hubbard parameter U has a similar effect 
as increasing the distance parameter in our analytical 
theory that means the increase of the Coulomb repulsin, 
see Eq. ([29)1 . Both effects push the centers of the wave 
functions from another. Notice, the parameter values 
used in the numerical simulations, correspond to rela- 
tively high non-adiabaticity of the system and strong an- 
harmonicity. Nevertheless, comparison of the three sets 
of figures corresponding to three different values of the 
Hubbard term in numerical simulations and Coulomb 
term in the analytical model shows that our analytical 
model gives rather good results even for rather strong 
electron repulsion. The approximations performed in de- 
riving our analytical model, in particular the transition 
to the continuum description and the expansions in de- 
riving the nonlinear potential and the separation of the 
center of mass motion from the relative motion, seem not 
to be decisive for the qualitative behavior of the bisolec- 
tron solution in the parameter range considered here. 



ative distance between the electrons becomes compara- 
ble with the lattice parameter for strongly localized elec- 
trons. Then, the singular Coulomb interaction (ex l/|a;|) 
should be replaced by the matrix elements with respect 
to the tight-binding orbitals and higher order derivatives 
in the continuum approximation should be taken into ac- 
count. 

The anharmonic bisoliton, found here and de- 
noted bisolectron (a soliton binding two electrons with 
Coulomb repulsion and Pauli's principle incorporated), 
can move with velocity up to the sound velocity, with 
both energy and momentum, maintaining finite values 
also at the sound velocity. Low-dimensional systems, 
mentioned in the Introduction, are characterized by the 
parameter values, for which the adiabatic approximation 
is valid. Their ground electron state is soliton-like and is 
extended over a few lattice sites. Therefore, we expect 
that in these systems pairing of electrons takes place at 
enough level of doping and that such bisolectron states 
are stable with account of the Coulomb repulsion and at 
the velocities up to the velocity of sound. Typical val- 
ues for poly-peptide macromolecules have been consid- 
ered that allow for the formation of a bisolectron bound 
state. 

Comparison of the energy of such a bisolectron with the 
energy of the two independent solitons binding a single 
electron each (two separate solectrons) shows that there 
is the gain of energy even including the Coulomb repul- 
sion. For illustration we have estimated the bisolectron 
binding energy for parameter values typical for biolog- 
ical macromolecules. We found that the lattice anhar- 
monicity significantly enhances pairing of electrons. We 
discuss the wave function of relative motion of two elec- 
trons in detail and include triplet states which lead to 
finite values of the Coulomb energy. Noteworthy is that 
there is a rather good qualitative agreement of the an- 
alytical model for two electrons in a singlet bisolectron 
state with account of Coulomb repulsion, developed here 
and the numerical simulations for the Hubbard model. 



VII. CONCLUSIONS 
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